function y = PARETO_inv (x, p)

   ## usage:  y = PARETO_inv (x, p)
   ##
   ## generalized Pareto inverse

   x = normcdf(x) ;
   y = p(2)/p(3) * ((1-x).^(-p(3)) - 1) + p(1) ;
   y = max(0, y) ;

endfunction
